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Abstract. In this paper we continue investigations of the 
Parker-shearing instability performing numerical simula- 
tions of the magnetic flux tube dynamics in the thin flux- 
tube approximation. We show that evolution of flux tubes 
resulting from numerical simulations is very similar to that 
of linear solutions if the vertical displacements are smaller 
than the vertical scale height H of the galactic disc. If the 
vertical displacements are comparable to H, the vertical 
growth of perturbations is faster in the nonlinear range 
than in the linear one and we observe a rapid inflation of 
the flux tube at its top, which leads to a singularity in 
numerical simulations, if only the cosmic rays are taken 
into account. Then we perform simulations for the case of 
nonuniform external medium, which show that the dom- 
inating wavelength of the Parker instability is the same 
as the wavelength of modulations of external medium. As 
a consequence of this fact, in the case of dominating cos- 
mic ray pressure, the dynamo a effect related to these 
short wavelength modulations is much more efficient than 
that related to the linearly most unstable long wavelengths 
modes of the Parker instability. Under the influence of dif- 
ferential forces resulting from differential rotation and the 
density waves, the a-effect is essentially magnified in the 
spiral arms and diminished in the intcrarm regions, what 
confirms our previous results obtained in the linear ap- 
proximation. 

Key words: Magnetic fields - Instabilities - Galaxies: 
magnetic fields - spiral - ISM: kinematics and dynamics 
of 



1. Introduction 

In the preceding paper (Hanasz and Lesch 1997, hereafter 
Paper I) we investigated the idea that the galactic dy- 
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namo action is driven by the Parker-shearing instability 
(Parker 1992, Hanasz & Lesch 1993, hereafter HL'93). We 
use the terms "Parker instability" or "Parker-shearing in- 
stability" depending on the context (i.e. on the current 
relevance of the rotational shearing effects). We presented 
the theory of slender magnetic flux tubes, analogous to 
that proposed by Spruit (1981), Spruit & Van Ballegooi- 
jen (1982) and Moreno-Insertis (1986), adopted to the case 
of galactic magnetic flux tubes. Then, performing the lin- 
ear stability analysis we calculated the dynamo coefficients 
resulting from the Parker-shearing instability of magnetic 
flux tubes in the thin flux tube approximation. The calcu- 
lations allowed us to derive a series of results which made 
conspicuous the power of the Parker-shearing instability 
as the mechanism generating the dynamo action (see Sect. 
4. of Paper I). 

It is obvious, however, that the linear approximation 
does not allow to inquire those aspect of the problem 
which are related to the evolution of perturbations for 
large amplitudes compared to the vertical scale height H. 
The problem results from the fact that the perturbations 
grow exponentially in the linear approximation and after a 
finite time period some physical quantities reach unphys- 
ical values. Since the possibilities of the extension of the 
work in the frame of analytical methods are rather limited, 
we decided to continue our studies by means of numerical 
simulations. We are going to present numerical solutions 
of the set of nonlinear equations describing the dynamics 
of flux tubes derived in Sect. 2 of Paper I together with 
the related dynamo transport coefficients. The solutions 
will allow us to verify the results of Paper I in the nonlin- 
ear regime and make a point to an old question about the 
final state of the Parker instability in galactic discs. 

The problem of the final state of the Parker instability 
with and without cosmic rays has been discussed inten- 
sively in the past. Parker (1968, 1969) argues that cos- 
mic rays leave the Galaxy by rapidly inflating bubbles on 
the magnetic field lines. This picture has been denied by 
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Mouschovias (1974, 1975) who has shown that: 1) with- 
out cosmic rays the Parker instability attains final, curved 
equilibrium states and 2) the final equilibria are attained 
also in the presence of cosmic rays if only the number 
of cosmic ray particles within a particular flux tube is 
conserved during the flux tube motion. In this model the 
cosmic ray pressure is uniform along the magnetic field 
line and varies with the volume of the flux tube. These 
2-dimensional solutions of Mouschovias have been tested 
with respect to instability against the perturbations in 
the 3-rd (radial) dimension by Asseo et al. (1978, 1980), 
Lcchieze-Rey et al. (1980) and appeared to be unstable 
even without cosmic rays. 

Some approaches to the Parker instability treated the 
cosmic ray gas as a relativistic gas, whose particles prop- 
agate with infinite speed along the magnetic field lines 
(see the discussion by Shu (1974) on this topic and ref- 
erences therein). It appears however, that the cosmic ray 
gas streaming along the magnetic field lines induces Alfven 
waves which exert a friction and slow the effective trans- 
port of cosmic rays down to the Alfven speed. Parker 
(1992) discusses his concept of inflating magnetic lobes, 
assuming explicitly that the inflation rate is limited by 
the Alfvenic streaming spend of the cosmic ray gas. 

The assumption of Mouschovias on the conservation 
of the number of cosmic ray particles within the flux tube 
cannot be retained because of the high production rate of 
cosmic rays due to the supernova explosions. There are 
at least two observational facts which strongly support 
the Parkers point of view. 1) The measured cosmic ray 
lifetime 10 7 yr in galactic discs implies that the cosmic 
rays should disconnect from the disc together with the 
trapping magnetic field. The characteristic time scale of 
Parker instability ~ 10 7 yr is consistent with this point 
of view. 2) In some edge-on galaxies a vertical magnetic 
field is observed in radio continuum (Hummel et al. 1991; 
Golla and Hummel 1994) extending in galactic halos. One 
should mention that the Parker instability is invoked to 
explain the formation of vertical magnetic fields in the 
galactic halo, among the other models like galactic foun- 
tains (Kahn and Brett, 1993) and galactic winds (Spencer 
and Cram, 1992). We are going to point out however, that 
the Parker instability with the contribution of cosmic rays 
can serve as a fully sufficient explanation. 

In HL'93 and Paper I we assumed that flux tubes are 
anchored at heavy molecular clouds. We have proposed 
that the typical intercloud distance determines the char- 
acteristic wavelength of the Parker instability. From Fig. 5 
of Paper I it follows that the range of unstable wavelengths 
depends on the magnetic pressure (via a). In the limit of 
weak magnetic field (e.g. a = 0.001) the marginal stability 
point is far below the wavelength 100 pc, then the typical 
distance between clouds (~ 100 pc) corresponds to an un- 
stable wavelength. Analogously, we can expect that even 
with the present magnetic field strength an instantaneous 
and local growth of the cosmic ray pressure due to super- 



nova explosions can help to overcome the magnetic field 
tension and shift the marginal stability point below given 
intercloud distance. In such a case the dominating wave- 
length of Parker instability can be much smaller than the 
frequently invoked value of a few vertical scalehcights. 

In the present section we generalize the physical setup 
and introduce a physically new effect - the propagation 
of flux tubes in an external medium which is nonuniform. 
The nonuniformities of the external medium are assumed 
to be spatially periodic and oscillating in time. It is our 
aim to check what is the evolution of a single flux tube in 
the nonuniform environment. The local and temporal ex- 
ternal gas condensations are intended to imitate the clouds 
in the ambient medium. We assume that such clouds ap- 
pear and disappear along the length of the flux tube mod- 
ifying the external density only, but have not a dynamical 
influence on the flux tube. This seems to be quite difficult 
at present to treat consistently the flux tube - cloud, or 
flux tube - flux tube collisions. Such a local and temporal 
external density modulations influence the flux tube inter- 
nal mass distribution via the magneto-hydrostatic balance 
condition (Eqn. (30) of Paper I) . Variations of the exter- 
nal density distribution force a continuous mass redistri- 
bution inside the tube triggering the possible instabilities 
all over the flux tube lifetime, not only at the initial mo- 
ment. In this respect our model seems to be more similar 
to of realistic interstellar medium than to the highly ide- 
alized uniform case. We are going to show that even a 
small amplitude external density modulations determine 
the characteristic wavelengths of the Parker modes. This 
in turn essentially influence the magnitude of the a-effect. 

The effects of the galactic dynamics on the Parker in- 
stability has been discussed in Paper I in the linear ap- 
proximation. The differential forces due to the axisymmet- 
ric differential rotation and the density waves influence es- 
sentially the dynamics of flux tubes as well as the dynamo 
transport coefficients. We found in Paper I that due to 
the presence of the radial differential force the curves of 
growth rate split into two separate parts at the wavenum- 
ber k cr it which is determined by the balance of the dif- 
ferential force, Coriolis force and the radial component of 
magnetic tension. In the range k cr it < k < k mar g (the 
Parker range), where k marg is the ordinary marginal sta- 
bility point, the linear solutions are stationary growing, 
as in the case without shear and the a-effect is essen- 
tially magnified. In the range < k < k crit (the shearing 
range) the linear solutions behave like propagating waves 
and the a-effect is essentially diminished. We pointed out 
that the density waves can modify the differential force 
periodically with a relatively large amplitude, so that a 
perturbation with a given wavenumber placed in vicinity 
of k cr it can move from the Parker to shearing range and 
vice versa depending on the phase of the density wave. 
Since the Parker range corresponds to the spiral arms and 
the shearing range to the interarm regions we showed that 
the a-effect can switch on in arms and -off in the interarm 
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regions. We predict that the regular component should 
be stronger in the interarm regions than in the arms due 
to the simultaneous action of strong shear and the mag- 
netic reconnection. These properties are very consistent 
with the observed arms of the uniform magnetic field in 
NGC 6946 (Beck and Hoernes 1 996) placed in between the 
optical arms. We are going to verify the above results in 
the nonlinear range by means of the nonlinear numerical 
simulations. 

The plan of this paper is as follows: In Section 2. wc 
present the applied numerical method, Section 3 is de- 
voted to the comparison of the numerical solutions with 
linear analytical solutions of Paper I. In Section 4 we an- 
alyze the final states of the Parker instability. In Section 
5 we calculate the temporal evolution of the dynamo co- 
efficients. In Section 6 the physical conditions are gener- 
alized by introducing modulations in the density of exter- 
nal medium, in which the flux tube propagate. In Section 
7 we examine the influence of the effects of differential 
forces in the nonlinear regime on the flux tube evolution 
and the dynamo action. Conclusions and possible relations 
of our theoretical results to the observational results are 
contained in Section 8. 



2. The method 

In Paper I we have derived the system of MHD equations 
governing evolution of a slender magnetic flux tube in the 
thin flux tube approximation. Geometry of the thin flux 
tube and the equations of motion have been described in 
the Section 2. of Paper I. The system of partial differ- 
ential equations describing evolution of the flux tube is 
given by equations (35), (36), (37) of Paper I for 9 un- 
known functions of time and position along the flux tube, 
together with the condition of magnetohydrostatic equi- 
librium of the flux tube with its surrounding (34). The 
unknown functions are: the Lagrangian coordinates of an 
element of the flux tube x, y, z, the components of the 
tangent vector L x , L y , L z and the components of veloc- 
ity v x , v y and v z . We shall not rewrite these equations 
here for brevity. Since we are going to solve the system 
of partial differential equations by means of ID numerical 
simulations we shall concentrate only on aspects related 
directly to this method. 

We shall use the numerical algorithm PDECOL taken 
from the NETLIB library (algorithm 540R, Madsen & 
Sincovec 1992). The package PDECOL is based on the 
method of lines and uses a finite element collocation pro- 
cedure (with piecewise polynomials as the trial space) for 
the discretization of the spatial variable x. The collocation 
procedure reduces the partial differential equation system 
to a semi-discrete system which then depends only on the 
time variable t. The time integration is then accomplished 
by use of standard techniques. The package solves the gen- 



eral system of nonlinear partial differential equations of at 
most second order of the form 

dU / TT dU d 2 U\ 

with given boundary conditions and an initial condition. 
In our approach we shall apply the periodic boundary con- 
ditions and the initial conditions of two types. The first 
type will be given by an eigenmode resulting from the lin- 
ear stability analysis. Such an initial state will allow to 
compare the linear analytical and the nonlinear numerical 
solutions. Results obtained with this initial state will serve 
for tests of the numerical algorithm. As the second type 
of the initial state we shall use a random superposition 
of periodic perturbations of the initial velocity and the 
vanishing initial displacements. Initial states of this type 
will allow for studies of the flux tube dynamics in con- 
ditions which more realistically resemble the nonuniform 
interstellar medium. 

3. The comparison of analytical and numerical so- 
lutions 

With the aim to test the correctness of the numerical solu- 
tions we shall compare them with the analytical linear so- 
lutions of the Paper I. First, wc shall apply the initial state 
of a small amplitude, which is given by the linear solution 
with a wavenumber k and an amplitude Zi n a. Then wc 
shall observe the temporal evolution of the numerical and 
analytical solutions which both start from the same initial 
state. We expect that before the nonlinear effects become 
important, these two solutions should evolve accordingly. 
We are going to compare the linear and the nonlinear am- 
plitudes of displacements of the flux tube gas in the three 
spatial directions at a given moment to. The linear so- 
lution will be represented by Xi = \Xi\ exp(cjjto), Yl = 
|Yi| exp(wjto) and Z\ = \Zi\exp(u)ito). The numerical 
nonlinear solution will be represented by the appropriate 
quantities X n = max \AX(s, to)\, Y n = max \AY(s, to)\ 
and Z n = max|AZ(s,i )|, where AX(s,t ), AY(s,t ) 
and AZ(s, to) are the Lagrangian displacements of an cl- 
ement of the flux tube with respect to its equilibrium po- 
sition at the time to. The maxima are taken over a full 
range of variability of the length parameter s. An exem- 
plary result obtained with the angular velocity typical for 
our Galaxy = CIq, the Oort constant A = 0, the ratios 
of cosmic ray and magnetic pressures to the gas pressure 
a = [3 = 1, the aerodynamic drag coefficient Cd — and 
the wavelengths of the perturbation H — 1330 pc (rep- 
resenting the most unstable mode) is shown in Fig. 1. 
Throughout the paper the initial height of the azimuthal 
flux tube is Zq = 100 pc, which is comparable to a half of 
the typical vertical scalchcight (H ~ 200 pc). The other 
parameters, which are not mentioned explicitly arc the 
same as those assumed in Paper I for our Galaxy. 
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Fig. 1. The comparison of the nonlinear numerical solution 
(full line) to the linear analytical solution (dotted line). The 
amplitudes of displacements in three directions vs. time are 
plotted. The basic parameters are: Q = Qa, A. = 0, a = /3 = 1, 
A = 1330 pc and C D = 
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Fig. 2. The comparison of the nonlinear numerical solution 
(full line) to the linear analytical solution (dotted line). The 
amplitudes of displacements in three dimensions vs. time are 
plotted. Parameters like in Fig. 1 but Cd = 1 and r = 10 pc. 



We note that the solutions are in a very good agree- 
ment especially concerning the coordinates x and y. How- 
ever, a small discrepancy in the vertical direction is re- 
markable. The reason for this discrepancy is in the lin- 
ear approximation of expression for the magnetic tension. 
Since in the linear approximation the curvature vector 
K has only the linear component, then the coefficient 
B 2 /Airpi contributes only in the zeroth order. This means 
that we take in linear approximation Bq instead of B 2 , 
which is overestimated at the upward displaced parts of 
the flux tube, inflated due to the cosmic rays. Then the 
vertical component of the magnetic tension is diminished 
in the nonlinear regime with respect to the linear approxi- 
mation. A similar argumentation applies also to the buoy- 
ancy term. This term is underestimated in the linear ap- 
proximation. This explains why the top of the flux tube 
rises faster in numerical simulations than in the linear ap- 
proximation. Thus, we can say that the numerical and 
analytical linear solutions are mutually consistent for rel- 
atively small vertical amplitudes Z < H . This allows us 
to expect that the basic conclusions of Paper 1, essential 
for the dynamo theory remain invariable in the nonlinear 
regime for Z < H . We shall demonstrate in the next sec- 
tions of this paper that this expectation is correct. More- 
over, we shall be able to take into account some effects, 
which have not been tractable within the frame of the 
linear approximation. 

As an example we shall examine the effect of aerody- 
namic drag force. HL'93 proposed that the aerodynamic 
drag force plays an important role for the dynamo trans- 
port coefficients, especially for the turbulent diffusivity. It 
was argued in that paper that the aerodynamic drag force 
slows down the vertical motions of fiuxtubes and thereby 
diminishes the vertical diffusion of magnetic field. The dis- 
cussion of the effect of the drag force has been omitted in 
Paper I since we focused our attention on the linear ap- 
proximation and the resulting properties. Solving the flux 
tube equations by means of the nonlinear numerical simu- 
lations we are able to take the drag force into account. A 
corresponding solution with the parameter set as in Fig. 1 
but with the aerodynamic drag force determined by the 
aerodynamic drag coefficient Cd = 1 and the flux tube 
radius r = 10 pc, is presented in Fig. 2. 

It is apparent in Fig. 2 that with the given values of pa- 
rameters, the aerodynamic drag force influences motion of 
the flux tube only moderately. One should be aware how- 
ever that the assumed values of the drag coefficient and the 
flux tube radius are essentially model dependent. More- 
over, the assumed physics underlying the aerodynamic 
drag effect may appear to be too simple for the case of 
galactic flux tubes. Nevertheless, one can notice that ac- 
cording to our expectations the drag force slows down the 
components of velocity perpendicular to the flux tube axis 
retaining the longitudinal velocity almost unchanged. 
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4. The final state of Parker instability 

In this section we shall continue to examine the nonlin- 
ear evolution of modes, which are introduced to the nu- 
merical experiment as the linear initial states, as it was 
described in the previous section. The first point, which 
deserves a special attention, is related to the final states of 
our simulations. We note that the simulations brake down 
in a certain moment T end . Moreover, we are not able to 
extend the simulations past T en d by making the grid res- 
olution more tiny. Typically we divide the spatial domain 
- the full period of the perturbation into 200 grid cells. 
Even the division of the spatial domain into 5000 cells 
does not allow to extend simulations past T en d. The more 
frequent sampling of the numerical solutions allows to ap- 
proach T en d closer and closer. In the following Figs. 3. and 
4. we illustrate what happens. We plot all the most rele- 
vant physical variables: 3 components of the Lagrangian 
displacement of the flux tube, together with the 3 compo- 
nents of velocity and acceleration, the internal density pi, 
the external density p e , the strength of magnetic field B 
inside the tube, the mass per unit length £, the flux tube 
radius r and the ratio of the total speed to the Alfven 
speed H/wa as a function of position along the tube s. 
Fig. 3. and Fig. 4. show these quantities at T = 125.0 Myr 
and T = 125.2 Myr, respectively (T end < 125.3 Myr ). 

We can easily notice that the internal density Pi ap- 
proaches at the top of flux tube at the moment directly 
preceding T en d, while the external density still remains 
nonvanishing. The vertical acceleration due to the buoy- 
ancy force starts to grow without limit within a limited 
time period. The appropriate panels of Figs. 3 and 4. ap- 
parently show this effect. All the other quantities behave 
regularly up to the moment T en d- It is worthwhile how- 
ever, to have a look at the other panels too. 

First of all, we notice that the topmost part of the tube 
starts to inflate (the flux tube radius r grows) at the time 
preceding T en d- We shall argue that the reason for this 
inflation is the cosmic ray pressure. We have assumed: 1) 
that the cosmic ray pressure is constant along the mag- 
netic field lines and 2) that it is fixed at zq (as if the flux 
tube were connected to an infinite reservoir of cosmic rays 
at z ). The gas density and pressure of the unperturbed 
disc decay exponentially with the increasing height z, thus 
at a certain height z\ the total (gas + magnetic + cosmic 
ray) pressure equals to the cosmic ray pressure at the ini- 
tial flux tube height zq. The cosmic ray pressure inside 
the flux tube element displaced from z to z\ is the same 
as the one at Zo- Then, this is not possible to balance 
the internal cosmic ray pressure if the flux tube clement 
is displaced further above z\, where the surrounding to- 
tal pressure is smaller than the cosmic ray pressure alone 
inside the raised element of the flux tube. This pressure 
imbalance has to result in a rapid inflation of the rising 
flux tube at its top. The inflation together with the out- 
flow of the thermal gas from the top cause that the internal 
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Fig. 3. The snapshot of the most relevant physical variables 
vs. position along the flux tube s, short before the brake down 
of the simulation. The basic parameters are: f2 = Q.g, A = 0, 
a = /3 = 1, A = 1330 pc and Cd = The units are: pc, 
Myr, lH-atom/cm 3 and pG. The plotted quantities are: the 
components of displacement, velocity and acceleration, internal 
and external densities pi, p e , the internal magnetic field, mass 
per unit length £, radius of the flux tube r and the ratio of the 
modulus of velocity to the Alfven speed |i>|/«a- 
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Fig. 4. The same as in previous figure 0.2 Myr later, but still 
before the break down. 
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gas density rapidly approaches zero. Since the buoyancy 
force is still nonvanishing and the internal mass density 
becomes infinitely small, the vertical acceleration has to 
become infinitely large. We propose to call the effect de- 
scribed above as a top singularity for brevity. 

There are however weak points of the presented sce- 
nario. The magnetic flux is conserved, the magnetic field 
diminishes together with the internal gas density in such a 
way that the internal Alfven speed grows, but the vertical 
acceleration grows even faster, so that the vertical speed 
of the top starts to exceed the internal Alfven speed (see 
Figs. 3. and 4., bottom-right panel). As it was already 
mentioned (Hartquist and Morfil, 1986; Parker 1992) the 
bulk streaming speed of the cosmic ray gas is limited to the 
Alfven speed approximately, thus the assumption on the 
infinite streaming speed is not fulfilled in the final stages 
of our simulations. The limited streaming speed is difficult 
to incorporate quantitatively in the actual considerations, 
but one can expect that this effect should drastically limit 
the inflation rate as soon as the top speed starts to exceed 
the Alfven speed (the cosmic rays are not able to follow 
the top). Finally, the top speed should saturate at the 
value of the order of Alfven speed, which is large anyway. 

On the other hand, the thin flux tube approximation 
becomes invalid in the numerical simulations as the flux 
tube diameter becomes comparable to the vertical scale 
lengths, what happens very close to the top singularity. 
The top singularity causes that we are not able to trace the 
evolution as well as the properties relevant to the dynamo 
theory after the moment T en d, even if all other parts of the 
flux tube in spite of its top (and its surrounding) behave 
regularly. 

Concluding this section we can say that the top sin- 
gularity is an important limitation for any numerical in- 
vestigations of the Parker instability which does not take 
into account the limited propagation speed of cosmic rays 
along the magnetic field lines. 

5. The dynamo coefficients 

The dynamo coefficients based on our numerical simula- 
tions can be calculated with the method applied in Paper 
I (formulae (82) and (84)). As a first step we shall present 
the temporal evolution of the helicity ad, diffusivity rjd, 
The magnetic Reynolds numbers R a , i? w and the dynamo 
number D for the parameter sets the same as that used in 
Figs. 1 and 2. This will allow us to make a first comparison 
of these coefficients resulting from numerical simulations 
to that obtained in the linear approximation in Fig. 4 of 
Paper I. 

Fig. 5 representing the case without aerodynamic drag 
shows a close similarity of the behaviour of the dynamo co- 
efficients to the results obtained in the frame of linear ap- 
proximation. This similarity has been naturally expected 
since the evolution of the shape is very similar in both the 
cases of numerical simulations and the linear approxima- 
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and D (full line) are shown in the lower panel. The aerody- 
namic drag coefficient Co = and the other parameters are 
like in Fig. 1. 
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Fig. 6. The same as in Fig. 5 in except of the aerodynamic 
drag coefficient, which is Cd = 1. This case corresponds to 
Fig. 2. of this paper. 
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tion. One minor difference is apparent. Since the nonlinear 
growth is faster than the linear one for large amplitudes 
of the Parker instability, the growth of the transport co- 
efficients ad and r\d is faster in the nonlinear case as well. 

The time dependence of the dynamo coefficients results 
from specific initial conditions taken into considerations: 
we begin with initial states of a small amplitude as de- 
scribed in the Section 3. We trace only the evolution of 
a single unstable flux tube. Since all the physical quanti- 
ties involved are time dependent, the dynamo coefficients 
are time dependent as well, fn more elaborative approach 
averaging over a statistical ensemble of flux tubes with 
a full spectrum of wavelengths and amplitudes of pertur- 
bations should be taken into account. This would allow 
to obtain saturated dynamo coefficients in galactic discs 
without density waves. 

The case of the aerodynamic drag coefficient (To = 1 is 
presented in Fig. 6. We notice that the aerodynamic drag 
diminishes the dynamo coefficients ad and r\d with respect 
to the case Cd = 0, because the motion of the flux tube 
is slowed down. 

It is now worthwhile to check if one of the main con- 
clusions of Paper I, concerning the a-effect in the weak 
magnetic field limit, is still valid in the nonlinear range. 
We shall examine the temporal evolution of the ad over 
the (a, /3)-plane. The simulations have been performed for 
a greed 10 times 5 of parameter points, covering the range 
0.0 < (3 < 1.0 and 0.0 < a < 0.5. The choice of the wave- 
length maximizing the ad coefficient has been made in the 
same way as in Paper I. All the simulations ended at the 
top singularity, while the top of the flux tube has been dis- 
placed by the distance approximately equal to the vertical 
scale length H . At the same time the values of the ad were 
always close to 1, sometimes even larger. The results are 
shown in Fig. 7, where we plot ad over the (a, /3)-plane 
for T — 40 and 70 Myr. In all the simulations the initial 
vertical displacement amplitude is 10 pc. 

We note that the results are very similar to those 
obtained in frame of the linear approximation. We can 
conclude that the strongest dynamo effect computed by 
means of the numerical simulations as well as in the lin- 
ear approximation is reached in a timcscale of a few tens 
of Myr for weak magnetic fields and the cosmic ray pres- 
sure comparable to the present value in the Milky Way, 
p cr ~ 0.5 • 10~ 12 dyn cm~ 2 . 
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Fig. 7. The temporal evolution of the ad coefficient over the 
(a, /3)-plane for T — 40 Myr (upper panel) and T — 70 Myr 
(lower panel). Because of the top singularity we replaced the 
missing points with the value ad = 1 km/s for T = 70 Myr. 

detailed description the problem looks as follows. Let us 
assign 

\ a = A(max(a d )), (2) 

for the wavelength maximizing the ad coefficient in the 
linear approximation, and A max for a certain very large 
wavelength for which the growth rate is close to maximum. 
The respective growth rates for these wavelengths are 



6. Nonlinear dynamics of magnetic flux tubes in a 
nonuniform external medium 

We have shown in Paper I, that in the limit of weak mag- 
netic field, the most significant contribution to the helicity 
ad is due to the short wavelength modes. But these modes 
are more stable than the long wavelength modes which 
contribute to ad only marginally (see Fig. 5 of Paper I). 
This leads us to the question if the generation process 
of magnetic field can still dominate over loses. In a more 



u ia = w t (max(a<()), w jmax = max(wj). (3) 

If we deal with a mixture of linear modes of various wave- 
lengths which start to grow with comparable initial am- 
plitudes at T — 0, then the main contribution to the total 
ad comes from the mode of the wavelength X a . In the lin- 
ear approximation ad grows with time proportionally to 
exp(2aji Q i). On the other hand, the main contribution to 
the turbulent diffusivity comes from the mode of the wave- 
length A max and grows proportionally to exp(2oJi max i) 
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with climax > L<Ji a . Then it is possible that in the case 
of a spectrum of unstable modes, with a variety of wave- 
length, the losses of magnetic field associated to the long 
wavelength modes will become dominant over the genera- 
tion process. 

To check this hypothesis in the nonlinear range, for 
the flux tube propagating in a uniform external medium 
we shall perform a numerical experiment. We assume that 
the initial state is a composition of periodic perturbations 
with wavenumbers given by 



kn — Tlkfl — — , 

Ao 



(4) 



where Ao is the maximum of wavelength determined by 
the size of the computational domain, and k n are the suc- 
cessive wavenumbers admitted by the periodic boundary 
conditions. The initial perturbations are constructed in 
such a way that only the velocity is perturbed at t = 0. 
The initial displacements are identically equal to 0. 



v(s, t = 0) = cos(fc n s + (f n ), 



(5) 



n=l 



where n max = 20 is the number of periodic components 
and ip n is a random phase given to each component for 
generality. We assume that each component has the same 
initial amplitude v n = 1 kms -1 and the initial phases are 
generated by the random number generator. We start sim- 
ulations with the above initial conditions and continue un- 
til the top singularity is reached. Then we make a plot of 
the shape of the fluxtube at the time directly preceding 
the top singularity. An example is shown in Fig. 8. We 
note that in spite of a compact region of strong deforma- 
tion due to the Parker instability, the rest of the flux tube 
remains relatively straight and close to the unperturbed 
equilibrium state. 

The temporal evolution of the dynamo coefficients for 
this case is presented in Fig. 9. We easily note differences 
with respect to the highly idealized case of the single com- 
ponent eigenmode initial perturbations. According to our 
expectations the dynamo coefficients a<2 and rjd are not 
proportional to the same common exponential factor in 
the case of uniform external medium. The turbulent dif- 
fusion coefficient grows in an exponential fashion appar- 
ently, while ad varies in a chaotic way around a value 
~ 0.01 kms -1 , which is rather small. For this reason R a 
is no longer close to a constant value, but diminishes with 
time. In consequence the dynamo number D diminishes 
much faster than R u . Then we note that the Parker in- 
stability would be very destructive for galactic magnetic 
field if the case discussed above was realistic. 

The above case is not realistic however, because the 
interstellar medium is very non-uniform. The presence of 
molecular clouds is the apparent signature of this fact. We 
expect that the clouds perturb permanently the motion of 
flux tubes. If we assume, for instance, that the flux tubes 
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Fig. 8. The shape of the flux tube perturbed with a super- 
position of periodic perturbations at the moment preceding 
appearance of the top singularity at T « 66 Myr 
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Fig. 9. The time dependence of the dynamo coefficients for the 
case presented in Fig. 8. Upper panel: ad (dashed line) and r/d 
(dotted line). Lower panel: R a (dashed line), R^ (dotted line) 
and D (full line). 
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are anchored to the molecular clouds as proposed by Beck 
et al. (1991), (see HL'93, Paper I) then the available wave- 
lengths of the Parker instability are related to the mean 
intcrcloud distance. The longer wavelengths of the order 
of galactic disc radius are suppressed due to the significant 
weight of molecular clouds and the escape of magnetic field 
by means of the very long wavelength modes is ineffective. 

Alternatively, one can assume that the flux tubes 
are infinitely long and uniform initially, but the ambient 
medium is cloudy. In the forthcoming considerations we 
shall apply a model of the flux tube - external medium 
interaction, resulting from the passage of the flux tube 
through a medium of a nonuniform, variable density. In 
order to limit a number of model assumptions we do not 
take into account the aerodynamic drag force, which in 
fact could very strongly influence the flux tube motion. We 
assume only that the density of external medium varies in 
a fashion resembling a standing wave, whose maxima ap- 
pear in variable places. Let us propose an analytical model 
of such variations taking 

Pe = Peo(z){l + ^(sm(k c s-0.5uj c t)-sm(k c s + 1.5Lu c t))(6) 

where the amplitude of the "cloud wave" A c = 0.2, and 
the wavelength A c = 2n/k c ~ 125 pc, equal to the most 
optimal wavelength for the ev-effect (see Fig. 5 of Paper 
I), which on the other hand is comparable to the typical 
intcrcloud distance 100 pc. The frequency of the external 
density modulations which we interpret as the flux tube - 
cloud collision frequency is 



Vc 



(7) 



where v c is the typical molecular cloud speed w 10 kms -1 . 
The above form of the external density modulations was 
incorporated in the subsequent simulation, where we ap- 
plied the same initial conditions as in the simulation lead- 
ing to the results presented in Figs. 8 and 9. The simu- 
lation had been continued again until the top singularity 
was reached at T ~ 39 Myr. The plot of the shape of 
fluxtube for this case is shown in Fig. 10. We notice that 
in the presence of the external density modulations, the 
whole length of the flux tube becomes buckled with the 
typical wavelength equal to the wavelength of the external 
density modulations, what is in contrast to the results ob- 
tained without the modulations. The temporal evolution 
of dynamo coefficients for this case is presented in Fig. 11. 

We note the essential differences with respect to the 
case without the modulations of external medium: 

1. The final value of the ad coefficient is significantly 
larger. Without external modulations we obtained 
ad ~ 0.05 kms , while with modulations it is ~ 
0.2 kms -1 . The final value of the % coefficient is ap- 
proximately the same as in the case without external 
modulations. 
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Fig. 10. The same as in Fig. 8, but with modulations of ex- 
ternal density. 
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Fig. 11. The time dependence of the dynamo coefficients for 
the case presented in Fig. 10. Upper panel: aa (dashed line) 
and r/d (dotted line). Lower panel: R a (dashed line), Ru, (dotted 
line) and D (full line). 
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2. After an initial transient period of about 20 Myr, the 
transport coefficients aa and rjd grow with a similar 
rate. This means that modes of the same growth rate 
contribute to the both a-effect and the turbulent dif- 
fusivity. 

3. The final value of the R a coefficient is « 0.5 while with- 
out modulations it was « 0.1. Let us recollect that R a 
contains the ratio ad/rjd, hence the generation rate of 
radial magnetic field is much more effective compared 
to losses in the presence of external modulations. 

We can summarize then the role of the modulations of ex- 
ternal medium as very favorable for the dynamo action. 
This is because the long wavelength modes, which are the 
most unstable in the linear approximation, (we consider 
now the weak magnetic field limit) do not dominate in the 
presence of the density modulations of external medium. 
The a-effect is strong due to the dominance of short wave- 
length modes. 

7. The effect of rotational shear 

In this section we present the results of numerical simu- 
lations taking into account the differential force resulting 
from the rotational shear (Fdi// = — 4AQp). As it has been 
described in Paper I the magnitude of shear is parameter- 
ized by the Oort constant A, which is equal to 0, -1/2 and 
-3/4 for rigid, flat and Keplcrian rotation curves respec- 
tively. The parameters of the simulations are: a = (3 = 1, 
f2 = f^c- We fix the wavenumber fco = fc C rit = 0.0044 pc" 1 
(A cr it = 1430 pc) for (A = —1/2). Whereas the wavenum- 
ber of perturbations is fixed, we shall vary the shearing 
parameter A around Aq = —1/2 corresponding to the flat 
rotation curve. In Paper I we have shown that the density 
waves introduce an additional differential force. The dif- 
ferential force corresponding to the linearity limit of the 
interstellar gas perturbations is equivalent to oscillations 
of A with the 'top-bottom' amplitude equal to about 50 % 
of |Ao| = 1/2. This linearity limit is related to a very small 
(1%) spiral density wave perturbation of the axisymmetric 
gravitational potential. 

Let us take for simplicity two fixed values A\ = —3/8 
and A 2 = —5/8 of the Oort constant A. The two val- 
ues are placed symmetrically around A and their differ- 
ence equals 50% of A$. For A\ the perturbation with fco 
is placed in the Parker range, while for A 2 in the shear- 
ing range. We shall perform two complementary runs of 
numerical simulations for these two values of the Oort con- 
stant. Our aim is to figure out the qualitative differences 
between these two cases, considering the basic physical 
variables already displayed in Figs. 3 and 4, and the dy- 
namo coefficients, as well. 

Let us focus on the shape of the flux tube and the 
related variables first. We shall analyze late stages of the 
evolution of the Parker-shearing instability, which precede 
the appearance of the top singularity by less than 5 Myr. 
The results are presented in Figs. 12 and 13. 
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Fig. 12. The snapshot of the most relevant physical vari- 
ables vs. position along the flux tube for the Parker range 
(A = Ai = -3/8, A = X crit = 1430 pc for A = 1/2). The 
other parameters and units are the same as in Fig. 3 
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In Fig. 12 (A = Ai, ko is in the Parker range) we note 
a regular, symmetric shape of the flux tube with apparent 
signatures of some nonlinear departures from the sinu- 
soidal shape. The central part of the azimuthal displace- 
ment X-curve is steepened, the top of Z-curve is vertically 
peaked and the bottom is shifted above the initial vertical 
position of the flux tube. In the linear regime the bottom 
displacement is negative and its absolute value is the same 
as the vertical displacement of the top. In contrast to the 
case without shear, the azimuthal velocity curve V x be- 
comes double periodic, what means that gas does not flow 
toward the lowest point of the flux tube at its lower parts, 
but instead it flows toward two points placed in between 
the top and the bottom. This is confirmed by the shape 
of the curve representing the mass per unit length of the 
flux tube £. We note two points of the gas concentration 
in between the positions of the top and bottom. The val- 
ues of pi approaching mean that the presented stage of 
the flux tube evolution is close to the top singularity and 
the assumptions of the thin flux tube approximation start 
to be violated near the top of the flux tube. The 
exceeding 1 means that the transport of cosmic rays along 
the magnetic field lines starts to be inefficient to inflate 
freely the topmost parts of the flux tube. 

In Fig. 13 (A = A 2l fc is in the shearing range) we 
note that the shape of the flux tube becomes irregular. 
We have demonstrated in Paper I that the solutions are 
the propagating waves in this range, so that the azimuthal 
position of the top is shifted with respect to its initial 
position at s = 0. The phase shift between the y- and 
z-displacements is less favorable for the magnitude of the 
ay coefficient than in the Parker range. In the linear range 
the formula for ad contains the quantity Im(Yi/Zi) (the 
imaginary part of radial to vertical displacement ampli- 
tude ratio), which is the largest for the phase difference 
about ±7r/2 between the vertical and horizontal displace- 
ments. Since this case represents the shearing range, ad is 
significantly reduced with respect to the Parker range as 
it is shown in Fig. 8 of Paper I. From the graph presenting 
the distribution of mass per unit length £ we learn that 
gas concentrates around a single point, which is maximally 
displaced in the positive y (radial) direction. Simultane- 
ously this is the lowest point on the flux tube. The point 
maximally displaced in the negative y direction has the 
lowest linear mass density and is the highest point of the 
flux tube. Similarly to the previously discussed case of 
A = A\, the presented stage of the flux tube evolution is 
at the limit of applicability of the used approximations. 

Comparing the above two cases we can say that the 
final mass distribution along the flux tube is a result of 
compromise between two different redistributing forces: 
the gravitational force, which tends to place gas in the low- 
est point of the flux tube and the differential force, which 
tends to place gas in the points which are extremely dis- 
placed in the radial direction. In the Parker range (the dif- 
ferential force is relatively weak), the two mass concentra- 



tions appear without a modification of the phase relations 
between the displacements in three spatial directions. For 
the stronger differential force in the shearing range the 
phase relations are modified and the two categories of con- 
centration points merge (i.e. the bottom becomes max- 
imally displaced in the radial direction). This effect di- 
minishes the dynamo a-effect forcing the disadvantageous 
phase relation between the radial and azimuthal displace- 
ments. It appears, however, that a residual a-effect re- 
mains in the shearing range, but it is much weaker than 
the a-effect in the Parker range. 

The temporal evolution of the dynamo coefficients for 
both the cases of Parker-shearing instability is presented 
in the subsequent Figs. 14 and 15, which are complemen- 
tary to Figs. 12 and 13 respectively. We note that ad is 
apparently larger in the Parker range (Fig. 14) than in 
the shearing range (Fig. 15), what confirms the previous 
results obtained in the linear regime. It is worth noting, 
however, that the dominance of the a-effect in the Parker 
range diminishes slightly during the nonlinear evolution 
of the Parker-shearing instability. The diffusivity r]d be- 
haves very similarly in both the cases because the growth 
rates are very similar. This implies the similarity of the 
magnetic Reynolds number R u . The magnetic Reynolds 
number R a and in consequence the dynamo number D 
are larger in the Parker range than in the shearing range 
following the same behaviour of ad- 
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Fig. 14. The time dependence of the dynamo coefficients for 
the Parker range (A = A\). Upper panel: ad (dashed line) and 
rj d (dotted line). Lower panel: R a (dashed line), (dotted 
line) and D (full line). 
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We can conclude that the nonlinear numerical simula- 
tions well confirm our previous linear results of Paper I, 
with the exception that the difference of between the 
Parker and shearing ranges becomes less pronounced in 
the nonlinear range. 

8. Discussion and Conclusions 

In this paper we present the results of numerical simula- 
tions solving the nonlinear system of partial differential 
equations, derived in Paper I for the dynamics of mag- 
netic flux tubes in galactic discs. The obtained numerical 
solutions appeared to be very consistent with the linear 
analytical solutions derived in Paper I for amplitudes of 
vertical displacements of the flux tubes smaller than the 
vertical scale length H . This consistency between the nu- 
merical and the analytical results allows us to confirm the 
previous conclusion: the weak magnetic field enables the 
fast growth of Parker instability and the strong dynamo 
action due to the weakness of magnetic tension if only the 
cosmic ray pressure is high enough. 

An essential advantage of the performed simulations 
follows from the fact that we have been able to trace the 
nonlinear evolution of the Parker unstable flux tubes. We 
have shown that due to the contribution of cosmic rays 
the top of the flux tube inflates rapidly at finite verti- 
cal displacements of the order of vertical scalehcight H. 
This gives rise to an enormous growth of the vertical buoy- 



ancy acceleration. This process is able to form the vertical 
structure of magnetic field lines extending to the galactic 
halo, as it is observed in some cases of edge-on galaxies 
like NGC 4631 (Hummel et al. 1991; Golla and Hummel 
1994) and M82 without any additional driving force. We 
note however that the inflation rate should be reduced if 
only the streaming speed of cosmic rays along the mag- 
netic field lines is limited, so that the velocities of vertical 
rise should saturate at the level comparable to the Alfven 
speed, typically several tens to hundreds of km/s. The 
last estimation depends on the particular values of the 
magnetic field strengths and gas densities, so that large 
variations of the Alfven speed are possible depending on 
galaxy. 

Reuter et al. (1992, 1994) point out that the halo 
of M82 is characterized by a filamentary structure with 
prominent gaps localized between these filaments. These 
features exhibit a preferential orientation normal to the 
plane of M82 and extend out to z-distanccs of ~ 1 kpc. 
The linear polarization up to 35 % was found there with 
the polarization degree increasing with z-distance. More- 
over, the vertical dependence of the spectral index of syn- 
chrotron radiation was interpreted as an evidence, that 
the relativistic plasma is streaming into the halo along 
vertical magnetic field lines with velocities in excess of 
1000 kms -1 . The M82 galaxy is known for its intense star 
formation activity at the center, what implies the intense 
production of cosmic rays, which are so important for our 
model. We can say that even if the Parker instability is not 
easily separable from the other effects like galactic winds 
and galactic fountains, the presented observational picture 
of the mentioned galaxies is consistent with our theoreti- 
cal model. We would like to point out that we can explain 
the vertical magnetic field structure of M82, the streaming 
of cosmic ray gas with very large vertical velocity as well 
as the generation of magnetic field within a single model 
of galactic dynamo driven by the Parker instability. 

We have made a point in this paper that if the cosmic 
ray pressure significantly exceeds the magnetic pressure, 
the Parker instability is able to produce an effective dy- 
namo action. This is possible if 1) the flux tubes are an- 
chored in molecular clouds with mean separation compa- 
rable to the most optimal wavelength for the a-effect (see 
Paper I), or 2) the flux tubes propagate in the nonuniform 
external medium as it has been discussed in this paper. 
Taking into account an exemplary model of weak mod- 
ulations of the external medium density we have shown 
that the dominating wavelengths of the buckling of flux 
tubes is the same as the wavelength of modulations of 
the external medium. We attribute these modulations to 
the presence of molecular clouds with typical separation 
of 100 pc. With these modulations the dynamo a-effect 
is essentially stronger than without the modulations and 
strong in comparison to the vertical diffusion of magnetic 
field. Without these modulations the diffusion rate tends 
to dominate over the generation rate. 
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Performing numerical simulations of the nonlinear evo- 
lution of flux tubes we have examined the role of the 
galactic differential forces resulting from the axisymmet- 
ric differential rotation for two distinct values of the Oort 
constant A. These simulations allowed us to figure out 
the qualitative differences between the Parker and the 
shearing ranges introduced in Paper I. Wc have found 
that for given perturbation with wavenumber fco, in the 
Parker range (weak shear) of the Parker-shearing insta- 
bility, the qualitative properties are similar to the case 
without shear. There is, however, an effect resulting from 
the presence of the radial differential force - the gas starts 
to gather not at the lowest point of the flux tube, but 
at two points extremely displaced in the radial direction. 
In contrast, we note a drastic change of the general be- 
haviour of the instability in the shearing range. The per- 
turbations become propagating waves. The gas starts to 
gather around one point, which is the lowest one and si- 
multaneously the most displaced in the radial direction. 
This qualitative change implies the modification of the 
phase difference between the radial and the vertical dis- 
placements. This effect in turn is responsible for the drop 
of the dynamo a^ coefficient. In the nonlinear regime the 
a-effect is still much more effective in the Parker range 
than in the shearing range, but the jump of its magni- 
tude is smaller than that obtained in the linear regime. 
Following Paper I we associate the Parker range with the 
enhanced a-effect with the spiral arms. Similarly we asso- 
ciate the shearing range with the diminished a-effect with 
the interarm regions. We note that a more elaborated ap- 
proach should take into account the higher cosmic ray 
production rate in spiral arms, which should additionally 
magnify the contrast of the a-effect between arms and the 
interarm regions. 

Summarizing, we can say that all the numerical sim- 
ulations well confirm our former results of Paper I, and 
additionally highlight the new nonlinear aspects. In Paper 
I we proposed a new model of galactic dynamo which is a 
cyclic process involving the action of the Parker-shearing 
instability, galactic differential rotation, density waves and 
magnetic reconnection. In this model the a-effect of the 
auj dynamo is essentially dependent on the phase of the 
galactic density wave and is localized in spiral arms. The 
complementary shear (the w-effect) is enhanced in the in- 
terarm regions. The dissipation process necessary for the 
dissipation of the small scale nonuniform magnetic field 
components is accomplished by the magnetic reconnec- 
tion. Our dynamo model is a kind of a fast dynamo pro- 
posed by Parker (1992), which is free of problems pointed 
out by Kulsrud and Anderson (1992). 

Our approach, based on nonlinear dynamics of flux 
tubes, allowed us to compute the dynamo transport coef- 
ficients with the back reaction of the large scale magnetic 
field taken into account. We found that the presence of 
cosmic rays and the rotational differential force allow to 
overcome magnetic tension, so that the dynamo transport 
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coefficients do not saturate under the influence of mag- 
netic field of the strength below the observed equipartition 
value. Thus, our model avoids the basic problems of mean 
field dynamos as described by Vainshtcin and Cattaneo 
(1992). 
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